Tarea 5: Bayesian Inference Part II

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Fecha límite de entrega:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Primera Parte: Preguntas Teóricas
  5. Segunda Parte: Elaboración de Código

Objetivo

a la uuuuultima tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la ultima parte del curso, los cuales se enfocan principalmente en aplicar inferencia bayesiana para generar regresiones lineales y estudiar métodos de obtención de la posterior mas poderosos, como es MCMC. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Documentación:

Primera Parte: Preguntas Teóricas

A continuación, se presentaran diferentes preguntas que abordan las temáticas vistas en clases. Por favor responda cada una de estas preguntas de forma breve, no más de 4 o 5 lineas.

Pregunta 1:

Señale algunos beneficios (no mas de dos) que nos brinda realizar una regresión lineal bayesiana sobre una regresión ajustada por mínimos cuadrados.

1.- La regresión Bayesiana permite entregar información previa (prior) 2.- Dados los resultados (posterior) es posible analizar de mejor manera el grado de incertidubre de estos.

Pregunta 2:

A continuación se presenta un modelo de regresión lineal bayesiana:

\[y_i \sim Normal(\mu, \sigma)\] \[\mu = \beta_0 + \beta_1*x\] \[\beta_0 \sim Normal(10,2)\] \[\beta_1 \sim Normal(0,1)\] \[\sigma \sim Uniform(0,50)\]

En base al modelo presentado, responda las siguientes preguntas:

  • Describa que representa cada una de las lineas del modelo.

  • Señale los parámetros que conforman a la regresión lineal.

  • Que objetivo cumple en el modelo lineal tener una distribución \(Normal(0,1)\) en el parámetro \(\beta_1\).

  • Que propiedad de las regresiones lineales nos entrega \(sigma\).

Respuesta:

Describa que representa cada una de las lineas del modelo.

  • \(y_i\) corresponde a la likelihood, la cual tiene distribución gaussiana.
  • \(\mu\) corresponde a una variable que representa al modelo lineal.
  • \(\beta_0\), \(\beta_1\) y \(\sigma\) corresponden a distribuciones esperadas, es decir, priors. \(\beta_0\) y \(\beta_1\) son gaussianas y \(\sigma\) uniforme.

Señale los parámetros que conforman a la regresión lineal

Los parámteros \(\beta_0\) y \(\beta_1\) son los interceptos de la regresión lineal. Por su parte, x corresponde a la entrada del modelo.

Que objetivo cumple en el modelo lineal tener una distribución \(Normal(0,1)\) en el parámetro \(\beta_1\)

El parámetro \(\beta_1\) indica la relación entre \(\mu\) y x, es por esto que se le entrega un rango de valores negativos y positivos, pasando por el cero, ya que existe la posibilidad de que la relación sea nula.

Que propiedad de las regresiones lineales nos entrega \(sigma\).

\(\sigma\) corresponde a la desviación estandar de la linalización, es decir, de \(\mu\).

Pregunta 3:

Explique de forma resumida como funciona el algoritmo de Laplace approximation utilizado para la regresión lineal. Señale el porque existe la necesidad de aplicar este modelo y los supuestos que se realizan con este método.

Respuesta:

Este método es conveniente mayormente para posteriores unimodales y simetrias, es decir, con distribuciónes aproximables por funciones gaussianas.

Laplace aproximation asume que el posterior tiene una distribución multivariable gaussiana. Esto se puede expresar de la siguiente forma:

\[ f(θ_1,..., θ_m|d) = N( \overrightarrow{\mu,} Σ). \] > Donde \(\overrightarrow{\mu,}\) corresponde al vector de medias y \(Σ\) a la matriz de covarianza.

Para este método no es necesario calcular el posterior normalizado, es decir, no se debe calcular f(d), lo cual es generalmente muy costoso de computar.

Pregunta 4:

Determine si las siguientes afirmaciones son verdaderas o falsas. Justifique las falsas.

  • El algoritmo de metropolis hasting solo funciona si la probabilidad de saltar de B a A es la misma que de A a B.
  • El algoritmo de Gibbs funciona en cualquier caso.
  • El algoritmo de metropolis hasting y de Gibbs son el mismo algoritmo, pero este ultimo es una variante del primero.

Respuesta Aquí

[F] El algoritmo de metropolis hasting solo funciona si la probabilidad de saltar de B a A es la misma que de A a B

El algoritmo de Metropolis trabaja con distribuciones simetricas, donde saltar de B a A es la misma que de A a B. Por su parte, el algoritmo de metropolis hasting es más general y permite distribuciones asimetrias del posterior.

[V] El algoritmo de Gibbs funciona en cualquier caso.

El algoritmo de Gibbs, al ser una variante del algoritmo de Metropolis hasting, permite distribuciones asimetrias del posterior.

[V] El algoritmo de metropolis hasting y de Gibbs son el mismo algoritmo, pero este ultimo es una variante del primero.

El algoritmo de Gibbs corresponde a una variante del algoritmo de Metropolis hasting. En esta variante se recorre un parámtero a la vez, igualando la “proposal distribution” al posterior.

Pregunta 5:

El algoritmo de Gibbs es mas eficiente que el de metropolis hasting. ¿Como se logra este efecto? ¿Existe alguna limitante de esta estrategia?

La eficiencia viene de recorrer un parámtero a la vez, igualando la “proposal distribution” al posterior en el parámetro de turno. Con esto, la “proposal distribution” siempre es aceptada.

Este método tiene dos principales limitaciones. Primero, los casos donde no se quiete trabajar con priors conjugados. Segundo, lo casos donde se tiene un gran número de variables.

Pregunta 6:

De una ventaja y una desventaja del algoritmo HMC.

Como desventaja se tiene un mayor costo computacional, devido al trabajo con gradientes. Como ventaja, se tiene que es má eficiente al obtener respuestas, y como consecuencia requiere de menos datos para poder obtener buenos resultados.

Pregunta 7:

Nombre y explique dos propiedades que son deseables en una cadena de Markov.

Se desea que una cadena de Markov tenga estacionalidad, es decir, converger a la posterior y variar en su vecindad. Además, se desea que sea aleatoria (well mixing), es decir, que no haya correlación entre los datos, lo cual se representa como tendencias en el tiempo.

Pregunta 8:

Explique cómo cross-validation, criterios de información y regularización ayudan a mitigar o medir los problemas de underfitting y overfitting.

Respuesta Aquí

Pregunta 9:

Diseñe una DAG para un problema causal inventado por usted de al menos 5 nodos (puede basarse en el ejemplo de Eugene Charniak) usando Dagitty y considere que la DAG tenga al menos: una chain, un fork y un collider. Muestre con dagitty todas las independencias condicionales de su DAG. Explique las independencias usando las reglas de d-separación.

Respuesta Aquí

#install.packages("dagitty")
library("dagitty")

my_dag <- dagitty('dag {
bb="0,0,1,1"
Caries [pos="0.644,0.745"]
Edad [pos="0.645,0.430"]
Frecuencia_lavado_dientes [pos="0.747,0.551"]
Ida_al_dentista [pos="0.379,0.428"]
Ingesta_azucar [pos="0.558,0.549"]
Proteccion_fluor [pos="0.389,0.638"]
Edad -> Frecuencia_lavado_dientes
Edad -> Ingesta_azucar
Frecuencia_lavado_dientes -> Caries
Ida_al_dentista -> Proteccion_fluor
Ingesta_azucar -> Caries
Proteccion_fluor -> Caries
}

')

plot(my_dag)

Segunda Parte: Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.6.3
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.1     v dplyr   1.0.6
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1
## Warning: package 'tibble' was built under R version 3.6.3
## Warning: package 'tidyr' was built under R version 3.6.3
## Warning: package 'readr' was built under R version 3.6.3
## Warning: package 'purrr' was built under R version 3.6.3
## Warning: package 'dplyr' was built under R version 3.6.3
## Warning: package 'forcats' was built under R version 3.6.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(dplyr)
library(tidyr)

# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 3.6.3
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
# Análisis bayesiano
library("rethinking")
## Loading required package: rstan
## Warning: package 'rstan' was built under R version 3.6.3
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 3.6.3
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:tidyr':
## 
##     extract
## Loading required package: parallel
## rethinking (Version 2.01)
## 
## Attaching package: 'rethinking'
## The following object is masked from 'package:purrr':
## 
##     map
## The following object is masked from 'package:stats':
## 
##     rstudent

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')

Pregunta 1: Regresión Lineal Bayesiana

El objetivo de esta pregunta es introducirlo en la aplicación de una regresión bayesiana. Con esto, buscaremos entender como calcular una regresión bayesiana en base al “motor” aproximación de Laplace, revisando como se comporta la credibilidad de sus predicciones y como la regresión lineal puede llegar a mutar a aplicar una transformación en el vector \(x\). Para responder esta pregunta centre su desarrollo solo en las clases y las funciones que nos brinda la librería rethinking.

Unos expertos en alometría deciden realizar un estudio de las alturas de unos niños en un colegio, buscando generar un regresor lineal bayesiano capaz de predecir la altura en base al peso de los alumnos. Para realizar este trabajo recopilan los datos table_height.csv, quien posee observaciones fisiológicas de 192 alumnos.

Parte I

En conocimiento los datos recopilados por los expertos, le solicitan realizar la siguiente serie de tareas:

  • Defina un modelo de regresión bayesiana, definiendo sus propios priors. Comente cada una de sus asunciones y señale a través de ecuaciones el modelo. Para definir los priors puede ser interesante la información recopilada en el siguiente link: Priors
  • Ajuste el modelo lineal utilizando el método de Laplace approximation. Estudie a través del comando precis los resultados obtenidos y coméntelos.
  • Gráfique el MAP de regresión lineal obtenida, añadiendo al gráfico el intervalo del \(95\%\) que se tiene sobre la media y los valores predecidos de la altura. Comente los resultados obtenidos y señale si su modelo logra ser un buen predictor de los valores estudiados.

Parte II

En base a los resultados obtenidos, el experto que trabaja con usted le señala que las alturas se suelen modelas con pesos logarítmicos, por lo que le sugiere añadir un logaritmo natural en el vector \(x\) que compone su modelo lineal. Realice nuevamente la regresión utilizando un intervalo del \(95\%\) sobre la media y los valores predecidos de la altura. Comente los resultados obtenidos, señalando si el modelo logra ajustar mejor los valores.

Respuesta:

Parte I

  • Defina un modelo de regresión bayesiana
d <- read.csv("table_height.csv", header = TRUE)
d.hw <- d[ , c("height","weight")]
summary(d.hw)
##      height           weight      
##  Min.   : 53.98   Min.   : 4.252  
##  1st Qu.: 89.13   1st Qu.:11.708  
##  Median :111.12   Median :16.981  
##  Mean   :108.32   Mean   :18.414  
##  3rd Qu.:127.72   3rd Qu.:23.417  
##  Max.   :158.12   Max.   :44.736

A continuación se define el modelo. La likelihood se asume como una gaussiana. El promedio de alturas es obtenido en base a el peso de la persona ponderado y sumado por los interceptos. Ambos interceptos tienen un prior gaussiano. Por un lado \(\beta_{0}\) debiese ser similar al promedio de las alturas, para esto se utiliza un promedio aproximado de las alturas de niños en el estudio revisado. Dado que existen un rango grande de alturas se elige una desviación de 10 en el promedio.

Por otro lado, \(\beta_{1}\) corresponde a la relación de la altura con el peso, tomando como base \(\beta_{0}\). Para este valor se escoge una normal, ya que esta abarca las relaciones más probables que estas pueden tener.

Por último, Se escoge una desviación entre 0 y 40 centimetros, esto dado la alta variedad de alturas en estos rangos de edades.

model1<-alist( 
    height ~ dnorm( mu, sigma ),
    mu <- b0 + b1*weight,
    b0 ~ dnorm( 125 , 10 ) ,
    b1 ~ dnorm( 0 , 1) ,
    sigma ~ dunif( 0 , 40 )
)
  • Ajuste el modelo lineal utilizando el método de Laplace approximation

Se obtiene un valor \(\beta_{0}\), tal como se esperaba, obtiene un valor cercano al promedio de los datos. Por su parte, \(\beta_{1}\) obtiene un valor que muestra la relación entre el peso y la altura, que tal como se esperaba es positiva. Por último, \(\sigma\) tiene un valor amplio, esto tiene sentido debido a que este estudio se realiza en edades de crecimiento.

b.reg1 <- quap(model1, data=d.hw)
precis(b.reg1, prob=0.95)
  • Gráfique el MAP de regresión lineal obtenida

A continuación se muestra la regreción lineal obtenida. Los resultados no logran predecir correctamente los valores estudiados. Aun así, dada la linealización, es posible apreciar que se minimiza el error en relación con los datos.

# samples from the posterior
post1 <- extract.samples( b.reg1, n= 1e4 )

weight.seq1 <- seq( from=0 , to=50 , by=1 )
mu.link1 <- function(weight) post1$b0 + post1$b1*weight
mu1 <- sapply( weight.seq1 , mu.link1 )
mu.mean1 <- apply( mu1 , 2 , mean )
mu.HPDI1 <- apply( mu1 , 2 , HPDI , prob=0.95 )

height.weight1 <- function(weight)
rnorm(
n=nrow(post1) ,
mean=post1$b0 + post1$b1*weight ,
sd=post1$sigma )

sim.height1 <- sapply( weight.seq1 , height.weight1)
height.HPDI1 <- apply( sim.height1 , 2 , HPDI , prob=0.95 )

plot( height ~ weight , data=d.hw , col=col.alpha(rangi2,0.5) )
lines( weight.seq1 , mu.mean1 )
shade( mu.HPDI1 , weight.seq1 )
shade( height.HPDI1 , weight.seq1 )

Parte II

A continuación se repite el procedimiento anterior pero esta vez utiliando pesos logarítmicos. Con este cambio se puede apreciar como la regresión se ajusta adecuadamente a los datos, tanto en forma como en la minimización del error.

model2<-alist( 
    height ~ dnorm( mu, sigma ),
    mu <- b0 + b1*log(weight),
    b0 ~ dnorm( 125 , 10 ) ,
    b1 ~ dnorm( 0 , 10) ,
    sigma ~ dunif( 0 , 40 )
)
b.reg2 <- quap(model2, data=d.hw)
precis(b.reg2, prob=0.95)
# samples from the posterior
post2 <- extract.samples( b.reg2, n= 1e4 )

weight.seq2 <- seq( from=0 , to=50 , by=1 )
mu.link2 <- function(weight) post2$b0 + post2$b1*log(weight)
mu2 <- sapply( weight.seq2 , mu.link2 )
mu.mean2 <- apply( mu2 , 2 , mean )
mu.HPDI2 <- apply( mu2 , 2 , HPDI , prob=0.95 )

height.weight2 <- function(weight)
rnorm(
n=nrow(post2) ,
mean=post2$b0 + post2$b1*log(weight) ,
sd=post2$sigma )

sim.height2 <- sapply( weight.seq2 , height.weight2)
height.HPDI2 <- apply( sim.height2 , 2 , HPDI , prob=0.95 )

# plot raw data
plot( height ~ weight , data=d.hw , col=col.alpha(rangi2,0.5) )
# draw MAP line
lines( weight.seq2 , mu.mean2 )
# draw HPDI region for line
shade( mu.HPDI2 , weight.seq2 )
# draw HPDI region for simulated heights
shade( height.HPDI2 , weight.seq2 )

Pregunta 2: MCMC

El objetivo de esta pregunta es lograr samplear, mediante la técnica de MCMC, la distribución gamma.

En general la distribución gamma se denota por \(\Gamma(\alpha,\beta)\) donde \(\alpha\) y \(\beta\) son parámetros positivos, a \(\alpha\) se le suele llamar “shape” y a \(\beta\) rate La densidad no normalizada de una distribución gamma esta dada por:

\[ f(x\mid \alpha,\beta) = \begin{cases} x^{\alpha -1}e^{-\beta x} ~ &\text{ si } x > 0\\ 0 ~&\text{si } x \leq 0 \end{cases} \]

  • Implemente el algoritmo de metropolis hasting utilizando \(q(\theta^{(t)} \mid \theta^{(t-1)}) = \mathcal{N}(\theta^{t-1},1)\), importante su función tiene que poder recibir:
    • La condición inicial \(\theta_{0}\).
    • Cantidad de repeticiones.
    • \(\alpha\) y \(\beta\).
    Escriba el algoritmo sin utilizar implementenaciones de la distribución gamma en r.

De ahora en adelante considere \(\alpha = 5\) y \(\beta = \frac{1}{5}\).

  • Considere \(\theta_{0} = 1\), grafique el histograma que obtiene para distintas cantidad de repeticiones, entre sus pruebas tiene que tener al menos una que tenga del orden de \(10^5\) repeticiones ¿Como cambia la distribución en funcion de las repeticiones?
  • Considere distintos valores de \(\theta_{0}\) y una cantidad de repeticiones grande (del orden de \(10^5\)), grafique las distribuciones que obtenga comente sus resultados ¿es importante la condición inicial del algoritmo?.
  • Utilizando \(\theta_{0}\) y cantidad de repeticiones conveniente (de acuerdo a lo que obtuvo en las partes anteriores) compare con la distribución real. Obs: En esta parte puede utilizar la distribución gamma que tiene implementado r para comparar con lo que usted realizo.

Respuesta:

gamma_fun <- function(x, alpha, beta){
  
  if (x > 0){
    (x ^ (alpha - 1)) * exp(-1*beta * x)
    
  }
  else{
    0
  }
}


MCMC_g_fun <- function(repetitions, omega_i, alpha, beta){
  
  p_omega_t_minus_1 <- omega_i
  #p_omega_a <- omega_i
  historia <- c()
  for(t in 1:repetitions){
    
    p_omega_a <- rnorm(1, mean=p_omega_t_minus_1, sd=1)

    # Suppose we have r
    gf_a <- gamma_fun(p_omega_a,         alpha, beta) 
    gf_b <- gamma_fun(p_omega_t_minus_1, alpha, beta)
    
    r_gf <- gf_b / gf_a #gf_a / gf_b
    
    r  <- r_gf
    #print("r")
    #print(r)
    if (rbinom(1, 1, min(1, r)) == 1){
      p_omega_t <- p_omega_a
    }
    else{
      p_omega_t <- p_omega_t_minus_1
    }
    
    p_omega_a <- p_omega_t
    historia <- append(historia, p_omega_a)
    
  }
  
  historia
}

P2 a)

library("ggplot2")

history <- MCMC_g_fun(1e2, 1, 5, 1/5)
history <- data.frame("h_data" = history)
ggplot(history, aes(x=h_data)) + 
  geom_density()

history <- MCMC_g_fun(1e3, 1, 5, 1/5)
history <- data.frame("h_data" = history)
ggplot(history, aes(x=h_data)) + 
  geom_density()

history <- MCMC_g_fun(1e4, 1, 5, 1/5)
history <- data.frame("h_data" = history)
ggplot(history, aes(x=h_data)) + 
  geom_density()

history <- MCMC_g_fun(1e5, 1, 5, 1/5)
history <- data.frame("h_data" = history)
ggplot(history, aes(x=h_data)) + 
  geom_density()
history <- MCMC_g_fun(1e5, 1, 5, 1/5)
history <- data.frame("h_data" = history)
ggplot(history, aes(x=h_data)) + 
  geom_density()

history <- MCMC_g_fun(1e5, 0.1, 5, 1/5)
history <- data.frame("h_data" = history)
ggplot(history, aes(x=h_data)) + 
  geom_density()

history <- MCMC_g_fun(1e5, 0.3, 5, 1/5)
history <- data.frame("h_data" = history)
ggplot(history, aes(x=h_data)) + 
  geom_density()

history <- MCMC_g_fun(1e5, 0.5, 5, 1/5)
history <- data.frame("h_data" = history)
ggplot(history, aes(x=h_data)) + 
  geom_density()

Respuesta Aquí

x <- seq(0, 4, 0.1)
y <- sapply(x, gamma_fun, alpha=5, beta=1/5)
mydf <- data.frame(x = x, res = y)

plot(x, y)



history <- MCMC_g_fun(1e4, 1, 5, 1/5)
history <- data.frame("h_data" = history)
ggplot(history, aes(x=h_data)) + 
  geom_density()
 

A work by CC6104